
# Declaration of Used Library

# For dataframe processing (shp and csv)
import pandas as pd
import geopandas as gpd
import numpy as np
from scipy.stats import pearsonr
from libpysal.weights import Queen

# For mapping geographical features
import osmnx as ox
import seaborn as sns
import contextily as ctx
from pysal.lib import weights
from pysal.explore import esda
from splot.esda import moran_scatterplot, lisa_cluster, plot_local_autocorrelation

# For plotting
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
from shapely.geometry import Point
%matplotlib inline
import mapclassify

%matplotlib inline

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/numba/core/decorators.py:262: NumbaDeprecationWarning: numba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.
  warnings.warn(msg, NumbaDeprecationWarning)
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/quantecon/lss.py:20: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def simulate_linear_model(A, x0, v, ts_length):
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/spaghetti/network.py:40: FutureWarning: The next major release of pysal/spaghetti (2.0.0) will drop support for all ``libpysal.cg`` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ``libpysal.cg`` geometries should prepare for this simply by converting to ``shapely`` geometries.
  warnings.warn(dep_msg, FutureWarning, stacklevel=1)


# Read CSV file, at desired level
initial_data = pd.read_csv("./data/shrug-antyodaya-csv/antyodaya_pc11subdist.csv")
# The content of the CSV will be explained on Part 02

# Check the data at a glimpse using head
initial_data.head()


# Review the number of records (row) and columns
print("Check the Number of Rows and Columns")
num_rows, num_columns = initial_data.shape
print(f"Number of rows: {num_rows}")
print(f"Number of columns: {num_columns}")

Check the Number of Rows and Columns
Number of rows: 5480
Number of columns: 171


# Get the name of the column (column title)
column_names = list(initial_data.columns)
print("The name of the columns of the data:" )

for i, column_name in enumerate(column_names, start=1):
    print(f"{i}. {column_name}")

The name of the columns of the data:
1. pc11_state_id
2. pc11_district_id
3. pc11_subdistrict_id
4. total_population
5. male_population
6. female_population
7. total_hhd
8. total_hhd_engaged_in_farm_activi
9. total_hhd_engaged_in_non_farm_ac
10. is_govt_seed_centre_available
11. availability_of_watershed_dev_pr
12. availability_of_rain_harvest_sys
13. availability_of_food_storage_war
14. availability_of_farm_gate_proces
15. availability_of_custom_hiring_ce
16. total_cultivable_area_in_hac
17. net_sown_area_in_hac
18. net_sown_area_kharif_in_hac
19. net_sown_area_rabi_in_hac
20. net_sown_area_other_in_hac
21. is_soil_testing_centre_available
22. is_fertilizer_shop_available
23. no_of_farmers_using_drip_sprinkl
24. area_irrigated_in_hac
25. total_unirrigated_land_area_in_h
26. availability_of_milk_routes
27. availability_of_poultry_dev_proj
28. availability_of_goatary_dev_proj
29. availability_of_pigery_developme
30. is_veterinary_hospital_available
31. availability_of_fish_community_p
32. availability_of_fish_farming
33. availability_of_aquaculture_ext_
34. total_hhd_with_kuccha_wall_kucch
35. total_hhd_have_got_pmay_house
36. total_hhd_in_pmay_permanent_wait
37. total_hhd_got_benefit_under_stat
38. total_hhd_in_perm_waitlist_under
39. piped_water_fully_covered
40. is_village_connected_to_all_weat
41. public_transport
42. availability_of_railway_station
43. total_hhd_having_pmsbhgy_benefit
44. availability_of_elect_supply_to_
45. availability_of_solor_wind_energ
46. total_hhd_having_solor_wind_ener
47. availability_of_panchayat_bhawan
48. csc
49. public_info_board
50. is_common_pastures_available
51. total_hhd_availing_pmuy_benefits
52. availability_of_public_library
53. recreational_centre
54. is_bank_available
55. is_bank_buss_correspondent_with_
56. is_atm_available
57. total_hhd_availing_pmjdy_bank_ac
58. is_post_office_available
59. telephone_services
60. is_broadband_available
61. is_pds_available
62. total_hhd_having_bpl_cards
63. availability_of_primary_school
64. is_primary_school_with_electrici
65. is_primary_school_with_computer_
66. is_primary_school_with_playgroun
67. is_primary_school_have_drinking_
68. availability_of_mid_day_meal_sch
69. total_primary_school_students
70. total_primary_school_teachers
71. availability_of_middle_school
72. availability_of_high_school
73. availability_of_ssc_school
74. no_of_children_not_attending_sch
75. availability_of_govt_degree_coll
76. total_grd_and_pg_in_village
77. is_vocational_edu_centre_availab
78. total_trained_trainee_under_skil
79. availability_of_jan_aushadhi_ken
80. total_hhd_registered_under_pmjay
81. is_community_waste_disposal_syst
82. total_hhd_with_clean_energy
83. is_community_biogas_waste_recycl
84. is_aanganwadi_centre_available
85. is_early_childhood_edu_provided_
86. total_childs_aged_0_to_3_years
87. total_childs_aged_0_to_3_years_r
88. total_childs_aged_3_to_6_years_r
89. total_childs_aged_0_to_3_years_i
90. total_childs_categorized_non_stu
91. total_anemic_pregnant_women
92. total_anemic_adolescent_girls
93. total_underweight_child_age_unde
94. total_male_child_age_bw_0_6
95. total_female_child_age_bw_0_6
96. total_minority_children_getting_
97. total_minority_hh_provided_bank_
98. no_of_physically_challenged_recv
99. total_hhd_with_more_than_2_child
100. availability_of_mother_child_hea
101. total_hhd_availing_pension_under
102. total_shg
103. total_hhd_mobilized_into_shg
104. total_no_of_shg_promoted
105. total_hhd_mobilized_into_pg
106. total_shg_accessed_bank_loans
107. is_bee_farming
108. is_sericulture
109. is_handloom
110. is_handicrafts
111. availability_of_community_forest
112. availability_of_minor_forest_pro
113. total_hhd_source_of_minor_forest
114. availability_of_cottage_small_sc
115. total_hhd_engaged_cottage_small_
116. availability_of_adult_edu_centre
117. total_no_of_registered_children_
118. total_no_of_children_0_to_6_year
119. total_no_of_pregnant_women
120. total_no_of_pregnant_women_recei
121. total_no_of_lactating_mothers
122. total_no_of_lactating_mothers_re
123. total_no_of_women_delivered_babi
124. total_no_of_children_in_icds_cas
125. total_no_of_young_anemic_childre
126. total_no_of_newly_born_children
127. total_no_of_newly_born_underweig
128. total_hhd_not_having_sanitary_la
129. total_no_of_eligible_beneficiari
130. total_no_of_beneficiaries_receiv
131. gp_total_no_of_eligible_benefici
132. gp_total_no_of_beneficiaries_rec
133. gp_total_hhd_eligible_under_nfsa
134. gp_total_hhd_receiving_food_grai
135. total_no_of_farmers_registered_u
136. total_no_of_farmers_subscribed_a
137. total_no_of_farmers
138. total_no_of_farmers_received_ben
139. total_no_farmers_adopted_organic
140. total_no_of_farmers_add_fert_in_
141. total_no_of_elected_representati
142. total_no_of_elect_rep_oriented_u
143. total_no_of_elect_rep_undergone_
144. total_approved_labour_budget_for
145. total_expenditure_approved_under
146. total_area_covered_under_irrigat
147. total_hhd_having_piped_water_con
148. pc11_land_area
149. canal
150. surface_water
151. ground_water
152. other_irrigation
153. any_primary_sch_toilet
154. mandi
155. regular_market
156. weekly_haat
157. phc
158. chc
159. sub_centre
160. closed_drainage
161. open_covered_drainage
162. open_uncovered_drainage
163. open_kuchha_drainage
164. fpo
165. no_electricity
166. internal_pucca_road
167. livestock_ext_services
168. _mean_p_miss
169. _core_p_miss
170. _target_weight_share
171. _target_group_max_weight_share


# Here, we will calculate the number of unique values in subdistrict and district column
# However, we need to drop the NA value of total_population
# because in my opinion, total population could not be easily replaced by certain value, such as mean
# It makes more sense if we change the value to density x total area in each district or subdistrict
# This is high effort, while we just need to see the glimpse at initial time.

def calculate_num_unique_values(df, column):
    # Filter out rows where 'total_population' is missing
    filtered_df = df.dropna(subset=['total_population'])
    
    # Calculate the number of unique values for the specified column
    num_unique_values = filtered_df[column].nunique()
    
    return num_unique_values

# Calculate the number of unique values for each specified column
num_unique_values_district = calculate_num_unique_values(initial_data, 'pc11_district_id')
num_unique_values_subdistrict = calculate_num_unique_values(initial_data, 'pc11_subdistrict_id')

# Print the number of unique values for each column
print("Number of unique values in pc11_district_id:", num_unique_values_district)
print("Number of unique values in pc11_subdistrict_id:", num_unique_values_subdistrict)

Number of unique values in pc11_district_id: 613
Number of unique values in pc11_subdistrict_id: 5444


selected_columns = ['pc11_state_id','pc11_district_id','pc11_subdistrict_id',
                    'total_population','male_population','female_population','total_hhd',
                    'total_grd_and_pg_in_village','total_shg',
                    'total_no_of_lactating_mothers','total_no_of_lactating_mothers_re',
                    'total_childs_aged_0_to_3_years','total_childs_aged_0_to_3_years_i'
                    ]

# Select only the columns in the 'selected_columns' array
subset_df = initial_data[selected_columns]


#Basically doing "for loops" and check every index in each of columns, is there any NA?
for x in selected_columns:
    if subset_df[x].isna().any() == True:
        na_loc = subset_df.index[subset_df[x].isna()].tolist()
        print(f"{x} with total {len(na_loc)} NaN data ({round(len(na_loc)/len(subset_df[x])*100,ndigits=2)}%), at index {na_loc}")

total_population with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
male_population with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
female_population with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_hhd with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_grd_and_pg_in_village with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_shg with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_no_of_lactating_mothers with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_no_of_lactating_mothers_re with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_childs_aged_0_to_3_years with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]
total_childs_aged_0_to_3_years_i with total 36 NaN data (0.66%), at index [62, 288, 1526, 1560, 1561, 1580, 1793, 3540, 4131, 4149, 4167, 4175, 4274, 4307, 4313, 4322, 4325, 4326, 4327, 4328, 4329, 4339, 4354, 4356, 4543, 4680, 4704, 4719, 4720, 4800, 4844, 4859, 4883, 4929, 5115, 5413]


# Remove NA value of point 1 and 2 above
columns_to_remove_na = ['total_shg', 'total_population']

# Drop rows with NA values in those columns
subset_df = subset_df.dropna(subset=columns_to_remove_na)


columns_to_remove_zero = ['total_no_of_lactating_mothers']

# Drop rows where specified columns contain only 0
subset_df = subset_df.loc[~(subset_df[columns_to_remove_zero] == 0).all(axis=1)]

# DO the same for total child ages 0-3 yrs
columns_to_remove_zero = ['total_childs_aged_0_to_3_years']
subset_df = subset_df.loc[~(subset_df[columns_to_remove_zero] == 0).all(axis=1)]


# No 01
subset_df['percent_immunized_child'] = (subset_df['total_childs_aged_0_to_3_years_i'] / subset_df['total_childs_aged_0_to_3_years']).round(3) * 100

# No 02
subset_df['percent_grd'] = (subset_df['total_grd_and_pg_in_village'] / subset_df['total_population']).round(3) * 100

# No 03
subset_df['percent_shg'] = (subset_df['total_shg'] / subset_df['total_population']).round(5)

# No 04
subset_df['percent_icds_women'] = (subset_df['total_no_of_lactating_mothers_re'] / subset_df['total_no_of_lactating_mothers']).round(3) * 100


# Recheck if there is NA value again
column_where_na = []
row_where_na = []

# Select all columns inside the subset_df DataFrame
all_columns = subset_df.columns

# Line below is used by me to check in each of column, if there is NA. then if there is NA, put it in the new empty array above.
for col in all_columns:
    if subset_df[col].isna().any():
        na_loc = subset_df.index[subset_df[col].isna()].tolist()
        print(f"{col}: NaN at index {na_loc}")
        column_where_na.append(col)
        row_where_na.extend(na_loc)

if not column_where_na:
    print("No NaN values found. Great!")
else:
    print("There are NaN values. Please check.")
    print("Columns with NaN values:", column_where_na)
    print("Rows with NaN values:", row_where_na)

No NaN values found. Great!


# I defined which column that I will plot in histogram
columns_hist = [
    'total_grd_and_pg_in_village',
    'total_shg',
    'total_no_of_lactating_mothers_re',
    'total_childs_aged_0_to_3_years_i'
]

# Then, I set the number of columns and rows for the layout
num_cols = 2
num_rows = 2

# Set up subplots for each column
fig, axes = plt.subplots(nrows=num_rows, ncols=num_cols, figsize=(10, 6))
axes = axes.flatten()

# Create histograms in each of the columns
for i, column in enumerate(columns_hist):
    axes[i].hist(subset_df[column].dropna(), bins=20, color='blue', alpha=0.7)
    axes[i].set_xlabel('Value')
    axes[i].set_ylabel('Frequency')

    if column == 'total_grd_and_pg_in_village':
        axes[i].set_title('Number of Graduates and Postgrads')
        axes[i].set_xlim(0, 50000)  # Set x-limit from 0 to 10000

    if column == 'total_shg':
        axes[i].set_title('Number of SHGs')
        axes[i].set_xlim(0, 8000)  # Set x-limit from 0 to 10000

    if column == 'total_no_of_lactating_mothers_re':
        axes[i].set_title('Number of Lactating Mothers Receiving ICDS')
        axes[i].set_xlim(0, 12000)  # Set x-limit from 0 to 10000

    if column == 'total_childs_aged_0_to_3_years_i':
        axes[i].set_title('Number of Children Age 0-3 Years with Immunization')
        axes[i].set_xlim(0, 40000)  # Set x-limit from 0 to 10000    

# Add a single title for all subplots
fig.suptitle('Histogram of Socioeconomic Key Variables in India', fontsize=16)
plt.tight_layout()
plt.show()


# Declareing which columns that we want to see the correlation
columns_corr = [
    'total_childs_aged_0_to_3_years_i',
    'total_grd_and_pg_in_village',
    'total_shg',
    'total_no_of_lactating_mothers_re'
]

# Create a subset dataframe with only the specified columns as specified above
subset_corr_df = subset_df[columns_corr]

# Change the name of the variable so that it is understandable by HUMAN!
subset_corr_df.columns = ['No. of Childs (0-3) Immunized','Education Level', 'No. of SHGs', 'No. of Lact. Mothers with ICDS']

# Create a pair plot
pair_plot = sns.pairplot(subset_corr_df, kind='scatter', diag_kind='hist')

# Change the colour to black (scatter plot)
for i, j in zip(*plt.np.triu_indices_from(pair_plot.axes, 1)):
    sns.scatterplot(x=subset_corr_df.columns[j], y=subset_corr_df.columns[i], data=subset_corr_df, color='black', ax=pair_plot.axes[i, j])

# Add orange linear regression lines with confidence intervals
for i, j in zip(*plt.np.triu_indices_from(pair_plot.axes, 1)):
    sns.regplot(x=subset_corr_df.columns[j], y=subset_corr_df.columns[i], data=subset_corr_df, scatter=False, color='orange', ax=pair_plot.axes[i, j])

# Add title
pair_plot.fig.suptitle('Pair Plot of Key Variables regarding Immunization Level of Childs (0-3) in India', y=1.02, fontsize=16)

# Show the plot
plt.show()


# Create a heatmap of the correlation matrix
correlation_matrix = subset_corr_df.corr()
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f', linewidths=.5)

# Show the plot
plt.title('Heatmap of Socio-economic Variables related to Immunization of Childs (0-3) in India', fontsize=12)
plt.show()


# Define which variables we want to observe
columns_name = [
    'percent_immunized_child',
    'percent_grd',
    'percent_shg',
    'percent_icds_women',
]


# Focus on dependent vs all independent variables on H1-H3
subset_corr_df = subset_df[columns_name]
subset_corr_df.columns = ['No. of Childs (0-3) Immunized','Education Level', 'No. of SHGs', 'No. of Lact. Mothers with ICDS']

# Create a pair plot with 'percent_immunized_child' on the y-axis and others on the x-axis
pair_plot = sns.pairplot(subset_corr_df, y_vars=['No. of Childs (0-3) Immunized'], x_vars=subset_corr_df.columns[1:], kind='scatter', diag_kind='hist')

# Show the plot
pair_plot.fig.suptitle('Pair Plot of Key Variables regarding Immunization Level of Childs (0-3) in India', y=1.12, fontsize=12)
plt.show()


# Create a heatmap of the correlation matrix
correlation_matrix = subset_corr_df.corr()
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f', linewidths=.5)

# Show the plot
plt.title('Heatmap of Socio-economic Variables (per capita) related to Immunization of Childs (0-3) in India', fontsize=12)
plt.show()


# Calculate correlation matrix with p-values
corr_matrix = subset_corr_df.corr()
p_values = np.zeros_like(corr_matrix)

for i in range(corr_matrix.shape[0]):
    for j in range(corr_matrix.shape[1]):
        _, p_value = pearsonr(subset_corr_df.iloc[:, i], subset_corr_df.iloc[:, j])
        p_values[i, j] = p_value

# Create a list of p-values
p_values_list = []

for i in range(corr_matrix.shape[0]):
    for j in range(corr_matrix.shape[1]):
        p_values_list.append((subset_corr_df.columns[i], subset_corr_df.columns[j], p_values[i, j]))

# Print the list of p-values
for p_value_info in p_values_list:
    print(f"P-value between {p_value_info[0]} and {p_value_info[1]}: {p_value_info[2]:.3f}")

P-value between No. of Childs (0-3) Immunized and No. of Childs (0-3) Immunized: 0.000
P-value between No. of Childs (0-3) Immunized and Education Level: 0.000
P-value between No. of Childs (0-3) Immunized and No. of SHGs: 0.000
P-value between No. of Childs (0-3) Immunized and No. of Lact. Mothers with ICDS: 0.000
P-value between Education Level and No. of Childs (0-3) Immunized: 0.000
P-value between Education Level and Education Level: 0.000
P-value between Education Level and No. of SHGs: 0.000
P-value between Education Level and No. of Lact. Mothers with ICDS: 0.000
P-value between No. of SHGs and No. of Childs (0-3) Immunized: 0.000
P-value between No. of SHGs and Education Level: 0.000
P-value between No. of SHGs and No. of SHGs: 0.000
P-value between No. of SHGs and No. of Lact. Mothers with ICDS: 0.000
P-value between No. of Lact. Mothers with ICDS and No. of Childs (0-3) Immunized: 0.000
P-value between No. of Lact. Mothers with ICDS and Education Level: 0.000
P-value between No. of Lact. Mothers with ICDS and No. of SHGs: 0.000
P-value between No. of Lact. Mothers with ICDS and No. of Lact. Mothers with ICDS: 0.000


map_path = './data/shrug-pc11subdist-poly-gpkg/subdistrict.gpkg'
map = gpd.read_file(map_path)
gdf = map


# Merge Data of subset_dataframe (subset_df) with the geo data frame (gdf)
# But first, let's examine the type of the column
gdf.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 5969 entries, 0 to 5968
Data columns (total 5 columns):
 #   Column               Non-Null Count  Dtype   
---  ------               --------------  -----   
 0   pc11_state_id        5969 non-null   object  
 1   pc11_district_id     5969 non-null   object  
 2   pc11_subdistrict_id  5969 non-null   object  
 3   subdistrict_name     5966 non-null   object  
 4   geometry             5969 non-null   geometry
dtypes: geometry(1), object(4)
memory usage: 233.3+ KB


# and the subset_df
subset_df.info()
subset_df.columns

<class 'pandas.core.frame.DataFrame'>
Index: 5420 entries, 1 to 5479
Data columns (total 17 columns):
 #   Column                            Non-Null Count  Dtype  
---  ------                            --------------  -----  
 0   pc11_state_id                     5420 non-null   int64  
 1   pc11_district_id                  5420 non-null   int64  
 2   pc11_subdistrict_id               5420 non-null   int64  
 3   total_population                  5420 non-null   float64
 4   male_population                   5420 non-null   float64
 5   female_population                 5420 non-null   float64
 6   total_hhd                         5420 non-null   float64
 7   total_grd_and_pg_in_village       5420 non-null   float64
 8   total_shg                         5420 non-null   float64
 9   total_no_of_lactating_mothers     5420 non-null   float64
 10  total_no_of_lactating_mothers_re  5420 non-null   float64
 11  total_childs_aged_0_to_3_years    5420 non-null   float64
 12  total_childs_aged_0_to_3_years_i  5420 non-null   float64
 13  percent_immunized_child           5420 non-null   float64
 14  percent_grd                       5420 non-null   float64
 15  percent_shg                       5420 non-null   float64
 16  percent_icds_women                5420 non-null   float64
dtypes: float64(14), int64(3)
memory usage: 762.2 KB

Index(['pc11_state_id', 'pc11_district_id', 'pc11_subdistrict_id',
       'total_population', 'male_population', 'female_population', 'total_hhd',
       'total_grd_and_pg_in_village', 'total_shg',
       'total_no_of_lactating_mothers', 'total_no_of_lactating_mothers_re',
       'total_childs_aged_0_to_3_years', 'total_childs_aged_0_to_3_years_i',
       'percent_immunized_child', 'percent_grd', 'percent_shg',
       'percent_icds_women'],
      dtype='object')


gdf.pc11_subdistrict_id = gdf.pc11_subdistrict_id.astype(float)
subset_df.pc11_subdistrict_id = subset_df.pc11_subdistrict_id.astype(float)


# Name the merged variable as geo_pop
geo_pop = gdf.merge(subset_df, on='pc11_subdistrict_id')
geo_pop.head()


# Because pc11_state_id_x = pc11_state_id_y and pc11_state_id_x = pc11_state_id_y
# We will delete column of y
# and rename the column pc11_state_id_x = pc11_state_id and pc11_state_id_x = pc11_state_id

geo_pop = geo_pop.drop(['pc11_state_id_y', 'pc11_district_id_y'], axis=1)

# Rename columns
geo_pop = geo_pop.rename(columns={'pc11_state_id_x': 'pc11_state_id', 'pc11_district_id_x': 'pc11_district_id'})


# Check if the crs is lose. Turns out that the crs is not lose, it still uses EPSG 4326. Seems fine, right?? :D
geo_pop.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich


mapdist_path = './data/shrug-pc11dist-poly-gpkg/district.gpkg'
map_dist = gpd.read_file(mapdist_path)


columns_to_sum = ['total_population', 'male_population', 'female_population', 'total_hhd',
                  'total_grd_and_pg_in_village', 'total_shg',
                  'total_no_of_lactating_mothers', 'total_no_of_lactating_mothers_re',
                  'total_childs_aged_0_to_3_years', 'total_childs_aged_0_to_3_years_i']

# Create a pivot table
pivot_table_df = subset_df.pivot_table(values=columns_to_sum, index='pc11_district_id', aggfunc='sum')

pivot_table_df


# Convert to astype
map_dist.pc11_district_id = map_dist.pc11_district_id.astype(float)
geo_pop.pc11_district_id = geo_pop.pc11_district_id.astype(float)
pivot_table_df.index = pivot_table_df.index.astype(float)

# Next, merge using the indices
geo_pop_district = map_dist.merge(pivot_table_df, left_on='pc11_district_id', right_index=True)
geo_pop_district.head()


# No 01
geo_pop_district['percent_immunized_child'] = (geo_pop_district['total_childs_aged_0_to_3_years_i'] / geo_pop_district['total_childs_aged_0_to_3_years']).round(3) * 100

# No 02
geo_pop_district['percent_grd'] = (geo_pop_district['total_grd_and_pg_in_village'] / geo_pop_district['total_population']).round(3) * 100

# No 03
geo_pop_district['percent_shg'] = (geo_pop_district['total_shg'] / geo_pop_district['total_population']).round(5)

# No 04
geo_pop_district['percent_icds_women'] = (geo_pop_district['total_no_of_lactating_mothers_re'] / geo_pop_district['total_no_of_lactating_mothers']).round(3) * 100


geo_pop_district.pc11_district_id = map_dist.pc11_district_id.astype(float)
pivot_table_df.index = pivot_table_df.index.astype(float)

# Merge based on pc11_district_id
geo_pop = pd.merge(geo_pop, geo_pop_district[['pc11_district_id', 'district_name']], on='pc11_district_id', how='left')


# Setup figure and axis
f, ax = plt.subplots(1)
# Add layer of polygons on the axis
gdf.plot(ax=ax)
# Add figure title
f.suptitle('India Map')
# Display
plt.axis('equal')
# Add labels to axes
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.show()


# We can also add the map from CTX! To make it more beautiful and you can easily imagine it above the atlas map
# Load local Authority geometries using pc11_subdistrict_id code as index
map = gpd.read_file(map_path).set_index('pc11_subdistrict_id', drop=False)
ax = map.plot(figsize=(6, 6), alpha=0.3, color='red');
# Add background map, expressing target CRS so the basemap can be
# reprojected (warped)
ctx.add_basemap(ax, crs=map.crs)
ax.set_title('Map of India')

Text(0.5, 1.0, 'Map of India')


number_cut = 1000
highest_imun = geo_pop.sort_values('percent_immunized_child',ascending=False).head(number_cut)
highest_shgs = geo_pop.sort_values('percent_shg',ascending=False).head(number_cut)
highest_grd = geo_pop.sort_values('percent_grd',ascending=False).head(number_cut)
highest_ICDS = geo_pop.sort_values('percent_icds_women',ascending=False).head(number_cut)


f, ax = plt.subplots(1, figsize=(6, 6))

background = geo_pop.plot(facecolor='blue', linewidth=0.025, ax=ax, alpha=0.2)
highest_ICDS.plot(alpha=0.9, facecolor='green', linewidth=0, ax=ax)
highest_imun.plot(alpha=0.7, facecolor='red', linewidth=0, ax=ax)

ax.set_axis_off()
f.suptitle('1000 Subdistricts with Highest Number of Lactating Mother (ICDS) and Children with Immunization in India')

# Create custom legend entries
legend_elements = [
    Line2D([0], [0], color='blue', alpha=0.2, lw=2, label='All Areas'),
    Line2D([0], [0], marker='o', color='green', label='Highest % Lact. Mother w ICDS'),
    Line2D([0], [0], marker='o', color='red', label='Highest % Child w Immunization'),
]
ax.legend(handles=legend_elements)
plt.axis('equal')
plt.show()


num_division = 11 # Divide into how many section
classi = mapclassify.Quantiles(geo_pop['percent_immunized_child'], k=num_division)
classi

Quantiles

    Interval       Count
------------------------
[  0.00,  47.30] |   497
( 47.30,  57.12] |   494
( 57.12,  63.40] |   490
( 63.40,  68.20] |   511
( 68.20,  72.50] |   490
( 72.50,  76.20] |   498
( 76.20,  80.10] |   492
( 80.10,  83.70] |   494
( 83.70,  87.60] |   498
( 87.60,  92.50] |   491
( 92.50, 100.00] |   492


# Set up the figure to explain the Quantiles that I made
f, ax = plt.subplots(1, figsize=(8, 5))
# Plot the kernel density estimation (KDE)
sns.kdeplot(geo_pop['percent_immunized_child'], shade=True)
# Add a blue tick for every value at the bottom of the plot (rugs)
sns.rugplot(geo_pop['percent_immunized_child'], alpha=0.5)
# Loop over each break point and plot a vertical red line
for cut in classi.bins:
    plt.axvline(cut, color='red', linewidth=0.75)
# Display image
ax.set_title('KDE for Immunized Children (0-3) Variable')    
plt.show()

/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/960233916.py:4: FutureWarning: 

`shade` is now deprecated in favor of `fill`; setting `fill=True`.
This will become an error in seaborn v0.14.0; please update your code.

  sns.kdeplot(geo_pop['percent_immunized_child'], shade=True)


# Create Choropleths using the 'cividis' colormap
geo_pop.plot(column='percent_immunized_child', scheme='QUANTILES', alpha=1, k=num_division, \
         cmap='cividis',  # Use cividis colormap
         edgecolor='w', linewidth=0.1, figsize=(8, 8))
         
# Add title with proper centering
plt.title('Choropleths of Immunized Children (0-3) in India, per Subdistricts', pad=20, loc='center')

# Add colorbar with proper positioning
cax = plt.axes([0.95, 0.1, 0.03, 0.8])  # [x, y, width, height]
sm = plt.cm.ScalarMappable(cmap='cividis', norm=plt.Normalize(vmin=geo_pop['percent_immunized_child'].min(), vmax=geo_pop['percent_immunized_child'].max()))
sm._A = []  # empty array for the data range
cbar = plt.colorbar(sm, cax=cax)

# Set colorbar label
cbar.set_label('Percentage of children aged 0-3 years immunized')

# Show the plot
plt.show()


# Here, I will show you Top 10
n_top_districts = 10

# Set the seaborn color palette to a colorblind-friendly palette
sns.set_palette("colorblind")

# Sort the DataFrame by 'percent_immunized_child' in descending order
sorted_df = geo_pop_district.sort_values(by='percent_immunized_child', ascending=False)

# Select the top districts based on the variable n_top_districts
top_districts = sorted_df.head(n_top_districts)

# Create the bar chart for 'percent_immunized_child'
fig, ax1 = plt.subplots(figsize=(12, 6))

color = sns.color_palette()[0]  # Use the first color from the palette
ax1.set_xlabel('District')
ax1.set_ylabel('Percent Immunized Child', color=color)
ax1.bar(top_districts['district_name'], top_districts['percent_immunized_child'], color=color, label='Percent Immunized Child')
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_ylim(0, 100)  # Set the y-axis limit for 'Percent Immunized Child'

# Create the line chart for 'percent_icds_women'
ax2 = ax1.twinx()
color = sns.color_palette()[1]  # Use the second color from the palette
ax2.set_ylabel('Percent ICDS Women', color=color)
ax2.plot(top_districts['district_name'], top_districts['percent_icds_women'], color=color, marker='o', linestyle='--', label='Percent ICDS Women')
ax2.tick_params(axis='y', labelcolor=color)
ax2.set_ylim(0, 100)  # Set the y-axis limit for 'Percent ICDS Women'

plt.title(f'{n_top_districts} Districts with Highest Percentage of Immunized Child and ICDS Women')
fig.tight_layout()
fig.legend(loc='upper left', bbox_to_anchor=(0.8, 0.2))  # Set the legend to a custom position
plt.show()


# Sort the DataFrame by 'total_population' in ascending order
geo_pop_district.sort_values(by='total_population').head(5)

# To ma
desired_column_order = ['pc11_state_id', 'pc11_district_id', 'district_name', 'geometry',
                         'female_population', 'male_population',
                         'total_childs_aged_0_to_3_years', 'total_childs_aged_0_to_3_years_i',
                         'total_grd_and_pg_in_village', 'total_hhd',
                         'total_no_of_lactating_mothers', 'total_no_of_lactating_mothers_re',
                         'total_shg', 'percent_immunized_child',
                         'percent_grd', 'percent_shg', 'percent_icds_women', 'total_population']

# Reorder the columns
geo_pop_district = geo_pop_district[desired_column_order]


selected_columns = ['pc11_state_id', 'pc11_district_id', 'district_name', 'geometry',
                     'total_population', 'percent_immunized_child']
geo_pop_district[selected_columns].sort_values(by='total_population').head()


# Filter rows where 'total_population' is greater than or equal to 10000
geo_pop_district = geo_pop_district[geo_pop_district['total_population'] >= 10000]


# Set up the figure and axes
plt.figure(figsize=(8, 4))

# Create a box plot
sns.boxplot(x='total_population', data=geo_pop_district, color='skyblue')

# Add title and labels
plt.title('Total Population by District')
plt.xlabel('Total Population')

# Show the plot
plt.show()


n_top_districts = 10

# Set the seaborn color palette to a colorblind-friendly palette
sns.set_palette("colorblind")

# Sort the DataFrame by 'percent_immunized_child' in descending order
sorted_df = geo_pop_district.sort_values(by='percent_immunized_child', ascending=False)

# Select the top districts based on the variable n_top_districts
top_districts = sorted_df.head(n_top_districts)

# Create the bar chart for 'percent_immunized_child'
fig, ax1 = plt.subplots(figsize=(12, 6))

color = sns.color_palette()[0]  # Use the first color from the palette
ax1.set_xlabel('District')
ax1.set_ylabel('Percent Immunized Child', color=color)
ax1.bar(top_districts['district_name'], top_districts['percent_immunized_child'], color=color, label='Percent Immunized Child')
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_ylim(50, 100)  # Set the y-axis limit for 'Percent Immunized Child'

# Create the line chart for 'percent_icds_women'
ax2 = ax1.twinx()
color = sns.color_palette()[1]  # Use the second color from the palette
ax2.set_ylabel('Percent ICDS Women', color=color)
ax2.plot(top_districts['district_name'], top_districts['percent_icds_women'], color=color, marker='o', linestyle='--', label='Percent ICDS Women')
ax2.tick_params(axis='y', labelcolor=color)
ax2.set_ylim(50, 100)  # Set the y-axis limit for 'Percent ICDS Women'

plt.title(f'{n_top_districts} Districts with Highest Percentage of Immunized Child and ICDS Women')
fig.tight_layout()
fig.legend(loc='upper left', bbox_to_anchor=(0.8, 0.2))  # Set the legend to a custom position
plt.show()


n_bottom_districts = 10

# Set the seaborn color palette to a colorblind-friendly palette
sns.set_palette("colorblind")

# Sort the DataFrame by 'percent_immunized_child' in ascending order to get the bottom districts
sorted_df_bottom = geo_pop_district.sort_values(by='percent_immunized_child', ascending=True)

# Select the bottom districts based on the variable n_bottom_districts
bottom_districts = sorted_df_bottom.head(n_bottom_districts)

# Create the bar chart for 'percent_immunized_child'
fig, ax1 = plt.subplots(figsize=(12, 6))

color = sns.color_palette()[0]  # Use the first color from the palette
ax1.set_xlabel('District')
ax1.set_ylabel('Percent Immunized Child', color=color)
ax1.bar(bottom_districts['district_name'], bottom_districts['percent_immunized_child'], color=color, label='Percent Immunized Child')
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_ylim(0, 100)  # Set the y-axis limit for 'Percent Immunized Child'

# Create the line chart for 'percent_icds_women'
ax2 = ax1.twinx()
color = sns.color_palette()[1]  # Use the second color from the palette
ax2.set_ylabel('Percent ICDS Women', color=color)
ax2.plot(bottom_districts['district_name'], bottom_districts['percent_icds_women'], color=color, marker='o', linestyle='--', label='Percent ICDS Women')
ax2.tick_params(axis='y', labelcolor=color)
ax2.set_ylim(0, 100)  # Set the y-axis limit for 'Percent ICDS Women'

plt.title(f'{n_bottom_districts} Districts with Lowest Percentage of Immunized Child and ICDS Women')
fig.tight_layout()
fig.legend(loc='upper left', bbox_to_anchor=(0.8, 0.2))  # Set the legend to a custom position
plt.show()


# Change 'percentage_immunized' to 'percent_immunized_child' in geo_pop dataframe
geo_pop = geo_pop.rename(columns={'percentage_immunized': 'percent_immunized_child'})

# Select a district with high percentage of immunization also with its neighbour
# To know the neighbour, I use maps here https://www.mapsofindia.com/districts-india/

selected_districts = ['Pathanamthitta', 'Alappuzha', 'Kottayam', 'Kollam', 'Idukki', 'Ernakulam', 'Thrissur', 'Theni', 'Tenkasi', 'Tirunelveli']
area1 = selected_districts  # I will put the above districts to different array, so in the end I can combine all selected districts!
filtered_geo_pop = geo_pop[geo_pop['district_name'].isin(selected_districts)]

# Input data from OSMX
# I use bbox that I input manually after knowing the boundary of the districts and neighbours
north, south, east, west = 10.4, 8.8, 77.4, 76.2
health_facil = ox.features_from_bbox(north, south, east, west, tags={"building": "hospital", "amenity": ["hospital", "clinic","doctors","social_facility"]})
#water = ox.features_from_bbox(north, south, east, west, tags={'natural': 'water'})

# Convert the features into GeoDataFrame
health_facil_gdf = gpd.GeoDataFrame.from_features(health_facil)

num_division = 20  # Divide into how many sections
classi = mapclassify.Quantiles(filtered_geo_pop['percent_immunized_child'], k=num_division)

fig, ax = plt.subplots(figsize=(8, 8))

filtered_geo_pop.plot(ax=ax, column='percent_immunized_child', scheme='QUANTILES', alpha=1, k=num_division, \
                      cmap='viridis',  # Use viridis colormap
                      edgecolor='w', linewidth=0.1, figsize=(8, 8))

# Add colorbar
cax = plt.axes([0.95, 0.1, 0.03, 0.8])  # [x, y, width, height]
sm = plt.cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=geo_pop['percent_immunized_child'].min(), vmax=geo_pop['percent_immunized_child'].max()))
sm._A = []  # empty array for the data range
cbar = plt.colorbar(sm, cax=cax)

# Set colorbar label
cbar.set_label('Percentage of Children (0-3 years) with Immunization')

health_facil_gdf.plot(ax=ax, color='red', markersize=10, label='Hospitals and Health Facility')  # Use dark red color
#water_gdf.plot(ax=ax, color='blue')


# Add legend
ax.legend()
ax.set_ylim(south-0.2, north+0.2)
ax.set_xlim(west-0.2, east+0.2)

ctx.add_basemap(ax, crs=geo_pop.crs)

ax.set_title('Immunization Percentage in Pathanamthitta Districts and surroundings')
# Show the plot
plt.show()

/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/mapclassify/classifiers.py:1592: UserWarning: Not enough unique values in array to form 20 classes. Setting k to 17.
  self.bins = quantile(y, k=k)
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/mapclassify/classifiers.py:1592: UserWarning: Not enough unique values in array to form 20 classes. Setting k to 17.
  self.bins = quantile(y, k=k)
/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/2884155885.py:43: UserWarning: Legend does not support handles for PatchCollection instances.
See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler
  ax.legend()


selected_districts = ['Anjaw', 'Lohit', 'Changlang', 'Namsai']
area2 = selected_districts  # I will put the above districts into a different array so that I can combine all selected districts!
filtered_geo_pop = geo_pop[geo_pop['district_name'].isin(selected_districts)]

north, south, east, west = 28.4, 27, 97.5, 95.5
health_facil = ox.features_from_bbox(north, south, east, west, tags={"building": "hospital", "amenity": ["hospital", "clinic", "doctors", "social_facility"]})

# Convert the features into GeoDataFrame
health_facil_gdf = gpd.GeoDataFrame.from_features(health_facil)

num_division = 20  # Divide into how many sections
classi = mapclassify.Quantiles(filtered_geo_pop['percent_immunized_child'], k=num_division)

fig, ax = plt.subplots(figsize=(10, 10))

filtered_geo_pop.plot(ax=ax, column='percent_immunized_child', scheme='QUANTILES', alpha=1, k=num_division, \
                      cmap='viridis',  # Use viridis colormap
                      edgecolor='w', linewidth=0.1, figsize=(8, 8))

# Add colorbar
cax = plt.axes([0.95, 0.1, 0.03, 0.8])  # [x, y, width, height]
sm = plt.cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=geo_pop['percent_immunized_child'].min(), vmax=geo_pop['percent_immunized_child'].max()))
sm._A = []  # empty array for the data range
cbar = plt.colorbar(sm, cax=cax)

# Set colorbar label
cbar.set_label('Percentage of Children (0-3 years) with Immunization')

health_facil_gdf.plot(ax=ax, color='red', markersize=10, label='Hospitals and Health Facility')  # Use dark red color

# Add legend
ax.legend()

ctx.add_basemap(ax, crs=geo_pop.crs)

ax.set_title('Immunization Percentage in Anjaw Districts and surroundings')
# Show the plot
plt.show()

/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/mapclassify/classifiers.py:1592: UserWarning: Not enough unique values in array to form 20 classes. Setting k to 17.
  self.bins = quantile(y, k=k)
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/mapclassify/classifiers.py:1592: UserWarning: Not enough unique values in array to form 20 classes. Setting k to 17.
  self.bins = quantile(y, k=k)


selected_districts = ['Bhandara', 'Nagpur', 'Gondiya', 'Seoni', 'Balaghat', 'Chhindwara', 'Garhchirolo', 'Chandrapur', 'Wardha']
area3 = selected_districts  # I will put the above districts into a different array so that I can combine all selected districts!
filtered_geo_pop = geo_pop[geo_pop['district_name'].isin(selected_districts)]

north, south, east, west = 22.7, 19.5, 81, 78.5
health_facil = ox.features_from_bbox(north, south, east, west, tags={"building": "hospital", "amenity": ["hospital", "clinic", "doctors", "social_facility"]})

# Convert the features into GeoDataFrame
health_facil_gdf = gpd.GeoDataFrame.from_features(health_facil)

num_division = 20  # Divide into how many sections
classi = mapclassify.Quantiles(filtered_geo_pop['percent_immunized_child'], k=num_division)

fig, ax = plt.subplots(figsize=(10, 10))

filtered_geo_pop.plot(ax=ax, column='percent_immunized_child', scheme='QUANTILES', alpha=1, k=num_division, \
                      cmap='viridis',  # Use viridis colormap
                      edgecolor='w', linewidth=0.1, figsize=(8, 8))

# Add colorbar
cax = plt.axes([0.95, 0.1, 0.03, 0.8])  # [x, y, width, height]
sm = plt.cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=geo_pop['percent_immunized_child'].min(), vmax=geo_pop['percent_immunized_child'].max()))
sm._A = []  # empty array for the data range
cbar = plt.colorbar(sm, cax=cax)

# Set colorbar label
cbar.set_label('Percentage of Children (0-3 years) with Immunization')

health_facil_gdf.plot(ax=ax, color='red', markersize=10, label='Hospitals and Health Facility')  # Use dark red color

# Add legend
ax.legend()
ctx.add_basemap(ax, crs=geo_pop.crs)

# Show the plot
ax.set_title('Immunization Percentage in Bhandara Districts and surroundings')
plt.show()

/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/203090821.py:32: UserWarning: Legend does not support handles for PatchCollection instances.
See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler
  ax.legend()


# Join each hospital+health facility with the corresponding subdistrict
# Here, I use "inner" (use intersection of the geometry) and "within" (choose that is inside) 
# The data that will be joined are: health facility gdf and the main dataframe (filtered based on districts)
selected_districts = area1 + area2 + area3 
filtered_geo_pop = geo_pop[geo_pop['district_name'].isin(selected_districts)]

joined_data = gpd.sjoin(health_facil_gdf, filtered_geo_pop, how="inner", op="within")

# Calculate the number of hospitals per area
hospital_counts = joined_data.groupby('subdistrict_name').size().reset_index(name='health_facil_count')

# Calculate the area of each subdistrict
calc_area = filtered_geo_pop[['subdistrict_name','district_name','geometry','percent_immunized_child']].copy()
calc_area = calc_area.to_crs('EPSG:32633')
calc_area['total_area'] = calc_area.geometry.area / 10**6

# Merge hospital_counts and district_areas to create the new array
health_facility_density_df = pd.merge(hospital_counts, calc_area[['district_name','subdistrict_name','percent_immunized_child','total_area','geometry']], on='subdistrict_name', how='outer')
health_facility_density_df = health_facility_density_df.fillna(0)
health_facility_density_df['health_facility/sq.km'] = health_facility_density_df['health_facil_count'] / health_facility_density_df['total_area']

column_order = ['district_name', 'subdistrict_name','health_facil_count', 'total_area','health_facility/sq.km','percent_immunized_child','geometry']
health_facility_density_df = health_facility_density_df[column_order]
health_facility_density_df.head()

/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3488: FutureWarning: The `op` parameter is deprecated and will be removed in a future release. Please use the `predicate` parameter instead.
  if await self.run_code(code, result, async_=asy):
/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/4145092681.py:7: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: None
Right CRS: EPSG:4326

  joined_data = gpd.sjoin(health_facil_gdf, filtered_geo_pop, how="inner", op="within")


from scipy.stats import pearsonr

# Calculate correlation coefficient and p-value
correlation_coefficient, p_value = pearsonr(health_facility_density_df['health_facility/sq.km'], health_facility_density_df['percent_immunized_child'])

# Print results
print(f'Correlation coefficient: {correlation_coefficient}')
print(f'P-value: {p_value}')

Correlation coefficient: -0.045525608575021206
P-value: 0.5737801079884187


# Scatter plot with regression line
plt.figure(figsize=(10, 6))
sns.regplot(x='health_facility/sq.km', y='percent_immunized_child', data=health_facility_density_df)
plt.title('Health Facility Density vs. Percentage of Immunized Children')
plt.xlabel('Health Facility Density (per sq. km)')
plt.ylabel('Percentage Children (0-3) Immunized')
plt.xlim(0, 0.001)  # Set x-axis limits
plt.grid(True)
plt.show()


variables = 'pc11_district_id'
weight_df = geo_pop_district

# Check for duplicate entries
if weight_df[variables].duplicated().any():
    # If duplicates exist, drop them or handle them as per your data requirements
    weight_df = weight_df.drop_duplicates(subset=variables)

# Set the index
weight_df = weight_df.set_index(variables, drop=False)

# Create spatial weights matrix
w_queen = Queen.from_dataframe(weight_df, idVariable=variables)

/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/3861510448.py:13: FutureWarning: `idVariable` is deprecated and will be removed in future. Use `ids` instead.
  w_queen = Queen.from_dataframe(weight_df, idVariable=variables)
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/libpysal/weights/weights.py:224: UserWarning: The weights matrix is not fully connected: 
 There are 25 disconnected components.
 There are 5 islands with ids: 71.0, 88.0, 528.0, 639.0, 640.0.
  warnings.warn(message)


w_queen.islands

[71.0, 88.0, 528.0, 639.0, 640.0]


weight_df.set_index('pc11_district_id')
weight_df = weight_df.drop(w_queen.islands)


w_queen = Queen.from_dataframe(weight_df, idVariable=variables)
w_queen.transform = 'R'

/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/3237309050.py:1: FutureWarning: `idVariable` is deprecated and will be removed in future. Use `ids` instead.
  w_queen = Queen.from_dataframe(weight_df, idVariable=variables)
/Users/mochammadmiftahulfahmi/anaconda3/envs/install_gds_stack/lib/python3.9/site-packages/libpysal/weights/weights.py:224: UserWarning: The weights matrix is not fully connected: 
 There are 20 disconnected components.
  warnings.warn(message)


queen_card = pd.Series(w_queen.cardinalities)
sns.distplot(queen_card, bins=10)
plt.title('Distribution of Cardinalities for Queen Contiguity')

/var/folders/h2/ws329tlx2tj_9vdbqdyj50v40000gn/T/ipykernel_13231/3785541745.py:2: UserWarning: 

`distplot` is a deprecated function and will be removed in seaborn v0.14.0.

Please adapt your code to use either `displot` (a figure-level function with
similar flexibility) or `histplot` (an axes-level function for histograms).

For a guide to updating your code to use the new functions, please see
https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751

  sns.distplot(queen_card, bins=10)

Text(0.5, 1.0, 'Distribution of Cardinalities for Queen Contiguity')


weight_df['w_percent_immunized_child'] = weights.lag_spatial(w_queen, weight_df['percent_immunized_child'])


# Check whether it works or not :D
weight_df[['pc11_district_id', 'district_name','percent_immunized_child', 'w_percent_immunized_child']].head(5)


w_queen.neighbors[468.0]

[477.0, 476.0, 469.0, 470.0]


weight_df['percent_immunized_child_std'] = (weight_df['percent_immunized_child'] - weight_df['percent_immunized_child'].mean()) / weight_df['percent_immunized_child'].std()


weight_df['w_percent_immunized_child_std'] = weights.lag_spatial(w_queen, weight_df['percent_immunized_child_std'])


# Setup the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values
sns.regplot(x='percent_immunized_child_std', y='w_percent_immunized_child_std', data=weight_df, ci=None)
# Add vertical and horizontal lines
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)
# Display
plt.title('Moran Plot of Standardized Childs Immunization Percentages')
plt.show()


mi = esda.Moran(weight_df['percent_immunized_child'], w_queen)
print("Moran I value = ")
mi.I

Moran I value =

0.586662400879172


moran_scatterplot(mi);


# Setup the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))
# Plot values
sns.regplot(x='percent_immunized_child_std', y='w_percent_immunized_child_std', data=weight_df, ci=None)
# Add vertical and horizontal lines
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)
plt.text(1.75, 0.5, "HH", fontsize=25)
plt.text(1.5, -1.5, "HL", fontsize=25)
plt.text(-2, 1, "LH", fontsize=25)
plt.text(-1.5, -2.5, "LL", fontsize=25)
# Display
plt.show()


lisa = esda.Moran_Local(weight_df['percent_immunized_child'], w_queen)


# Break observations into significant or not
weight_df['significant'] = lisa.p_sim < 0.05
# Store the quadrant they belong to
weight_df['quadrant'] = lisa.q


weight_df['significant'].head(20)

pc11_district_id
468.0    False
469.0    False
470.0    False
471.0    False
472.0    False
473.0    False
474.0    False
475.0    False
476.0    False
477.0    False
478.0    False
479.0    False
480.0    False
481.0    False
482.0    False
483.0    False
484.0    False
485.0    False
486.0    False
487.0     True
Name: significant, dtype: bool


weight_df['quadrant'].head()

pc11_district_id
468.0    3
469.0    3
470.0    3
471.0    1
472.0    1
Name: quadrant, dtype: int64


lisa_cluster(lisa, weight_df);


# In more detailed graph

# Setup the figure and axis
f, ax = plt.subplots(1, figsize=(6, 6))

# Plot insignificant clusters
ns = weight_df.loc[weight_df['significant']==False, 'geometry']
ns.plot(ax=ax, color='k')

# Plot HH clusters
hh = weight_df.loc[(weight_df['quadrant']==1) & (weight_df['significant']==True), 'geometry']
hh.plot(ax=ax, color=sns.color_palette("coolwarm", as_cmap=True)(0.8))  # Red

# Plot LL clusters
ll = weight_df.loc[(weight_df['quadrant']==3) & (weight_df['significant']==True), 'geometry']
ll.plot(ax=ax, color=sns.color_palette("coolwarm", as_cmap=True)(0.2))  # Blue

# Plot LH clusters
lh = weight_df.loc[(weight_df['quadrant']==2) & (weight_df['significant']==True), 'geometry']
lh.plot(ax=ax, color=sns.color_palette("coolwarm", as_cmap=True)(0.5))  # Greenish-Blue

# Plot HL clusters
hl = weight_df.loc[(weight_df['quadrant']==4) & (weight_df['significant']==True), 'geometry']
hl.plot(ax=ax, color=sns.color_palette("coolwarm", as_cmap=True)(0.35))  # Light Red

# Style and draw
f.suptitle('LISA for Immunization Percentage on Children 0-3 yrs in India', size=14)
f.set_facecolor('0.75')
ax.set_axis_off()

plt.show()


plot_local_autocorrelation(lisa, weight_df, 'percent_immunized_child');


if mi.p_sim < 0.05:
    print("Spatial autocorrelation is statistically significant.")
    if mi.I > 0:
        print("Positive spatial autocorrelation (HH or LL).")
    else:
        print("Negative spatial autocorrelation (LH or HL).")
else:
    print("No statistically significant spatial autocorrelation detected.")

Spatial autocorrelation is statistically significant.
Positive spatial autocorrelation (HH or LL).

	pc11_state_id_x	pc11_district_id_x	pc11_subdistrict_id	subdistrict_name	geometry	pc11_state_id_y	pc11_district_id_y	total_population	male_population	female_population	...	total_grd_and_pg_in_village	total_shg	total_no_of_lactating_mothers	total_no_of_lactating_mothers_re	total_childs_aged_0_to_3_years	total_childs_aged_0_to_3_years_i	percent_immunized_child	percent_grd	percent_shg	percent_icds_women
0	24	468	3722.0	Lakhpat	MULTIPOLYGON (((68.39398 23.44476, 68.39264 23...	24	468	62117.708950	32088.750366	30025.930963	...	394.599993	566.165208	227.071608	204.869050	3082.118619	2619.901747	85.0	0.6	0.00911	90.2
1	24	468	3723.0	Rapar	MULTIPOLYGON (((70.53440 23.46489, 70.53222 23...	24	468	201902.906291	105105.320376	96796.567185	...	2531.541403	339.236735	930.099517	861.844679	6842.802254	4941.854062	72.2	1.3	0.00168	92.7
2	24	468	3724.0	Bhachau	MULTIPOLYGON (((70.45008 23.01226, 70.44904 23...	24	468	176504.000000	92678.000000	83460.000000	...	2577.000000	234.000000	952.000000	790.000000	6897.000000	4373.000000	63.4	1.5	0.00133	83.0
3	24	468	3725.0	Anjar	POLYGON ((70.23631 23.40849, 70.23631 23.40545...	24	468	150105.112093	79896.579777	70092.396031	...	3643.775958	491.505351	1147.882751	1027.598741	5673.050161	4004.627988	70.6	2.4	0.00327	89.5
4	24	468	3726.0	Bhuj	POLYGON ((69.78433 23.99110, 69.79143 23.98750...	24	468	238834.277866	119211.778880	113669.738389	...	2819.277177	230.560895	629.484862	570.504168	7436.929336	5231.051377	70.3	1.2	0.00097	90.6

	female_population	male_population	total_childs_aged_0_to_3_years	total_childs_aged_0_to_3_years_i	total_grd_and_pg_in_village	total_hhd	total_no_of_lactating_mothers	total_no_of_lactating_mothers_re	total_population	total_shg
pc11_district_id
1	23523.138329	26510.879626	2761.074569	2590.694210	704.340235	8981.478930	505.056064	479.194760	5.004162e+04	1.521253
2	360076.838747	392652.138478	29465.549515	19313.974815	33116.408459	124419.391360	4322.707964	3710.070887	7.551146e+05	1291.782192
3	45598.474752	47615.448440	2837.299473	2422.499434	3636.436542	16684.780804	791.205308	756.308394	9.321936e+04	350.046346
4	78150.007679	78037.294461	7344.751470	5050.990314	4134.530494	23630.458275	1106.910964	961.424367	1.760563e+05	653.763285
5	230263.244779	240528.135141	18368.934333	11508.827636	11365.975014	92762.747369	5829.314660	4037.745077	4.717603e+05	1120.292371
...	...	...	...	...	...	...	...	...	...	...
632	450007.906654	447742.145346	18915.735027	15588.340247	40968.391785	289061.923978	4448.443895	4225.381399	9.028997e+05	4561.162494
633	596281.419949	595956.683533	23808.052542	19714.351280	67835.559509	391963.629490	4917.232443	4590.071292	1.195993e+06	4454.343346
638	3154.866472	3307.155473	257.055194	213.852640	70.204149	1883.631334	22.681341	22.681341	6.462022e+03	35.642107
639	54314.627363	57265.113737	3669.293400	3063.362231	5526.490359	29831.195841	662.187411	636.354232	1.116362e+05	743.045473
640	53502.795148	58770.591779	3692.091661	3419.911097	4963.049693	26783.402756	739.984115	678.350943	1.123320e+05	725.685726

	pc11_state_id	pc11_district_id	district_name	geometry	female_population	male_population	total_childs_aged_0_to_3_years	total_childs_aged_0_to_3_years_i	total_grd_and_pg_in_village	total_hhd	total_no_of_lactating_mothers	total_no_of_lactating_mothers_re	total_population	total_shg
0	24	468.0	Kachchh	MULTIPOLYGON (((70.45008 23.01226, 70.44904 23...	6.918732e+05	7.659721e+05	52789.062307	36124.140136	29509.874064	332446.747616	7194.496554	6452.874545	1.465846e+06	5262.200653
1	24	469.0	Banas Kantha	MULTIPOLYGON (((71.24964 24.20926, 71.24207 24...	1.490841e+06	1.603543e+06	137018.185689	93865.516320	103408.588403	607873.050328	19340.773617	16878.197941	3.273774e+06	13546.601715
2	24	470.0	Patan	MULTIPOLYGON (((71.42507 23.96967, 71.42497 23...	5.204153e+05	5.651527e+05	49457.172519	30594.736482	64222.702638	258319.762048	9340.013570	7797.987739	1.089937e+06	14351.448633
3	24	471.0	Mahesana	POLYGON ((72.79975 24.07615, 72.80022 24.07529...	7.639744e+05	8.327545e+05	59425.341078	49149.072387	123971.352675	380288.478687	8296.862992	6994.736144	1.617866e+06	15483.493324
4	24	472.0	Sabar Kantha	POLYGON ((73.14784 24.47759, 73.14773 24.47410...	1.096367e+06	1.183334e+06	87560.815945	66651.606962	152636.838216	582154.513255	16705.055796	14374.734389	2.320364e+06	20658.175179

	pc11_state_id	pc11_district_id	district_name	geometry	total_population	percent_immunized_child
423	12	257.0	Dibang Valley	POLYGON ((96.11495 29.39999, 96.14150 29.36881...	292.000000	62.5
61	22	405.0	Janjgir - Champa	POLYGON ((82.64385 22.21103, 82.64583 22.20764...	2038.000000	100.0
62	22	406.0	Bilaspur	POLYGON ((82.08579 23.11360, 82.08697 23.11342...	2363.000000	10.3
637	35	638.0	Nicobars	MULTIPOLYGON (((93.71364 7.21016, 93.71297 7.2...	6462.021944	83.2
451	15	285.0	Serchhip	POLYGON ((92.89217 23.56077, 92.89312 23.55949...	15921.602747	96.7

	district_name	subdistrict_name	health_facil_count	total_area	health_facility/sq.km	percent_immunized_child	geometry
0	Chhindwara	Amarwara	4.0	2994.261732	0.001336	70.6	POLYGON ((8136551.045 4848193.458, 8136662.527...
1	Gondiya	Amgaon	1.0	1105.078742	0.000905	71.0	POLYGON ((8404800.713 4832755.934, 8404840.361...
2	Gondiya	Arjuni Morgaon	3.0	3585.795810	0.000837	82.4	POLYGON ((8447580.202 4719914.218, 8448069.824...
3	Balaghat	Baihar	7.0	9699.919428	0.000722	74.2	POLYGON ((8362487.172 5026853.348, 8362637.339...
4	Balaghat	Balaghat	17.0	4077.470261	0.004169	66.0	POLYGON ((8270767.000 4931971.672, 8271153.673...

The program will be structured as below¶

Part 01: Initial Data Check ¶

Part 02: Research Formulation ¶

Part 03: Data Cleaning ¶

Part 04: Spatial Visualization ¶

Part 05: Exploratory Spatial Data Analysis ¶

Part 06: Conclusion and Reflection ¶

	pc11_state_id	pc11_district_id	pc11_subdistrict_id	total_population	male_population	female_population	total_hhd	total_hhd_engaged_in_farm_activi	total_hhd_engaged_in_non_farm_ac	...	open_uncovered_drainage	open_kuchha_drainage	no_electricity	internal_pucca_road	livestock_ext_services	_target_weight_share	_target_group_max_weight_share
0	3	41	225	778.0	407.0	371.0	155.0	0.0	0.0	...	1.0	0.0	0.00000	0.00000	0.0	1.0	1.0
1	12	245	1563	446.0	216.0	230.0	95.0	40.0	45.0	...	0.0	1.0	0.00000	0.00000	0.0	1.0	1.0
2	12	250	1626	603.0	253.0	350.0	129.0	110.0	3.0	...	0.0	0.0	0.45937	0.54063	0.0	1.0	1.0
3	12	253	1679	1380.0	828.0	552.0	228.0	120.0	70.0	...	0.0	1.0	0.00000	0.00000	1.0	1.0	1.0
4	21	371	2774	297.0	155.0	142.0	70.0	55.0	15.0	...	0.0	0.0	0.00000	1.00000	0.0	1.0	1.0